# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "devtools", "sdmTMB", "sdmTMBextra", "terra", "mapplots",
"viridis", "visreg", "modelr", "future", "kableExtra")
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Import some plotting functions
# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")
# devtools::install_github("seananderson/ggsidekick") # not on CRAN
library(ggsidekick)
theme_set(theme_sleek())
# Set path
home <- here::here()
# For crossvalidation: paralell processing
plan(multisession)
set.seed(99) Fit stomach content models
Load libraries
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/main-analysis/03-fit-diet-models_cache/html"))Read old and new data, combine!
df <- read_csv(paste0(home, "/data/clean/stomachs.csv")) |>
mutate(depth_sc = (depth - mean(depth))/sd(depth),
year_f = as.factor(year),
month_f = as.factor(month),
ices_rect = as.factor(ices_rect),
pred_length_sc = (pred_length - mean(pred_length)) / sd(pred_length)) |>
filter(pred_length > 10) |>
filter(FR_spr < 0.4) |> # Check this by year...
filter(FR_sad < 0.4) |>
filter(FR_her < 0.4)
# Read the prediction grid...
pred_grid <- bind_rows(read_csv(paste0(home, "/data/clean/pred_grid_(1_2).csv")),
read_csv(paste0(home, "/data/clean/pred_grid_(2_2).csv")))
# Scale with respect to data!
df |> group_by(month) |> summarise(n = n()) |> arrange(desc(n))# A tibble: 8 × 2
month n
<dbl> <int>
1 3 4620
2 11 1842
3 12 1508
4 2 839
5 4 332
6 10 34
7 6 32
8 5 2
pred_grid <- pred_grid |>
drop_na(oxy, temp, sal, depth) |>
filter(quarter == 4) |> # Not needed in theory for saduria...
mutate(depth_sc = (depth - mean(df$depth)) / sd(df$depth),
year_f = as.factor(year),
month_f = as.factor(3),
year = as.integer(year),
pred_length_sc = 0,
ices_rect = as.factor(ices_rect)) |>
droplevels()Plot data!
df |>
pivot_longer(c(FR_spr, FR_her, FR_sad)) |>
ggplot(aes(year, value)) +
geom_jitter(height = 0, alpha = 0.5) +
coord_cartesian(ylim = c(0, 0.1)) +
facet_wrap(~name, ncol = 1) +
geom_smooth()Set up a mesh
mesh <- make_mesh(df, c("X", "Y"), cutoff = 6)
ggplot() +
inlabru::gg(mesh$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = df, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.3, vjust = 3) +
labs(x = "Easting (km)", y = "Northing (km)")Spatial and temporal cross validation, and AIC, to compare spatial vs non spatial models
5-fold spatial cross validation & AIC to compare two models: spatial fields or ices rectangle as random effects
# Set up spatial clusters
clust <- kmeans(df[, c("X", "Y")], 5)$cluster
# Plot
df$clust_id <- clust
plot_map +
geom_point(data = df, aes(X*1000, Y*1000, color = as.factor(clust_id))) +
scale_color_viridis(discrete = TRUE) +
geom_sf(size = 0.1)Warning: Removed 33 rows containing missing values (`geom_point()`).
Warning: colourbar guide needs continuous scales.
# Sprat
# ICES as a random effects
spr_space_1_cv <- sdmTMB_cv(
FR_spr ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
spatial = "off",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Spatial random field
spr_space_2_cv <- sdmTMB_cv(
FR_spr ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
spatial = "on",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Herring
# ICES as a random effects
her_space_1_cv <- sdmTMB_cv(
FR_her ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
spatial = "off",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Spatial random field
her_space_2_cv <- sdmTMB_cv(
FR_her ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
spatial = "on",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Saduria
# ICES as a random effects
sad_space_1_cv <- sdmTMB_cv(
FR_sad ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
spatial = "off",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Spatial random field
sad_space_2_cv <- sdmTMB_cv(
FR_sad ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
spatial = "on",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
10-fold temporal cross validation & AIC to compare two models: spatial fields or ices rectangle as random effects
# Set up temporal clusters
# This is not very straightforward! Since with temporal clusters I don't have much spatial overlap
# Note we remove year from the model because it can't be estimated for single years
# clust <- as.numeric(as.factor(df$year))
#
# # Plot
# df$clust_id <- clust
#
# df <- df |> mutate(clust_id = round(clust_id/2))
#
# sort(unique(df$clust_id))
#
# df <- df |> mutate(clust_id = ifelse(!year == 2021, 1, 2))
#
# plot_map_fc +
# geom_point(data = df, aes(X*1000, Y*1000, color = year)) +
# scale_color_viridis() +
# facet_wrap(~clust_id) +
# geom_sf(size = 0.1)
#
# # Sprat
# # ICES as a random effects
# spr_temp_1_cv <- sdmTMB_cv(
# FR_spr ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
# data = df,
# mesh = mesh,
# spatial = "off",
# family = tweedie(link = "log"),
# fold_ids = clust,
# k_folds = length(unique(clust))
# )
#
# # Spatial random field
# spr_temp_2_cv <- sdmTMB_cv(
# FR_spr ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f),
# data = df,
# mesh = mesh,
# spatial = "on",
# family = tweedie(link = "log"),
# fold_ids = clust,
# k_folds = length(unique(clust))
# )
#
#
# # Herring
# # ICES as a random effects
# her_temp_1_cv <- sdmTMB_cv(
# FR_her ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
# data = df,
# mesh = mesh,
# spatial = "off",
# family = tweedie(link = "log"),
# fold_ids = clust,
# k_folds = length(unique(clust))
# )
#
# # Spatial random field
# her_temp_2_cv <- sdmTMB_cv(
# FR_her ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f),
# data = df,
# mesh = mesh,
# spatial = "on",
# family = tweedie(link = "log"),
# fold_ids = clust,
# k_folds = length(unique(clust))
# )
#
#
# # Saduria
# # ICES as a random effects
# sad_temp_1_cv <- sdmTMB_cv(
# FR_sad ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
# data = df,
# mesh = mesh,
# spatial = "off",
# family = tweedie(link = "log"),
# fold_ids = clust,
# k_folds = length(unique(clust))
# )
#
# # Spatial random field
# sad_temp_2_cv <- sdmTMB_cv(
# FR_sad ~ depth_sc + s(pred_length_sc, k=3) + (1|month_f),
# data = df,
# mesh = mesh,
# spatial = "on",
# family = tweedie(link = "log"),
# fold_ids = clust,
# k_folds = length(unique(clust))
# )AIC
# ICES as a random effects
fit_spr_m1 <- sdmTMB(
FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "off",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_spr_m1)✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
fit_her_m1 <- sdmTMB(
FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "off",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_her_m1)✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
fit_sad_m1 <- sdmTMB(
FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "off",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_sad_m1)✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
# Spatial random field
fit_spr_m2 <- sdmTMB(
FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_spr_m2)✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
fit_her_m2 <- sdmTMB(
FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_her_m2)✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✖ `ln_tau_G` standard error may be large
ℹ `ln_tau_G` is an internal parameter affecting `sigma_G`
ℹ `sigma_G` is the random intercept standard deviation
ℹ Try simplifying the model, adjusting the mesh, or adding priors
✖ `sigma_G` is smaller than 0.01
ℹ Consider omitting this part of the model
✔ Range parameter doesn't look unreasonably large
fit_sad_m2 <- sdmTMB(
FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
Print AIC & cross validation table
# Higher is better!
# likelihood from spatial cv
sum_loglik_space <- data.frame(Prey = c("Sprat", "Herring", "Saduria"),
rec = c(spr_space_1_cv$sum_loglik / 10000000,
her_space_1_cv$sum_loglik / 10000000,
sad_space_1_cv$sum_loglik / 10000000),
space = c(spr_space_2_cv$sum_loglik / 10000000,
her_space_2_cv$sum_loglik / 10000000,
sad_space_2_cv$sum_loglik / 10000000)) |>
mutate(rec_temp = rec, space_temp = space,
rec = ifelse(rec_temp > space_temp, paste0("**", rec, "**"), rec),
space = ifelse(space_temp > rec_temp, paste0("**", space, "**"), space)) |>
rename("sum loglik<sub>ICES rect</sub>" = rec,
"sum loglik<sub>spatial</sub>" = space) |>
dplyr::select(-rec_temp, -space_temp)
# Higher is better!
# likelihood from temporal cv
# sum_loglik_temp <- data.frame(Prey = c("Sprat", "Herring", "Saduria"),
# rec = c(spr_temp_1_cv$sum_loglik / 10000,
# her_temp_1_cv$sum_loglik / 10000,
# sad_temp_1_cv$sum_loglik / 10000),
# space = c(spr_temp_2_cv$sum_loglik / 10000,
# her_temp_2_cv$sum_loglik / 10000,
# sad_temp_2_cv$sum_loglik / 10000)) |>
# mutate(rec_temp = rec, space_temp = space,
# rec = ifelse(rec_temp > space_temp, paste0("**", rec, "**"), rec),
# space = ifelse(space_temp > rec_temp, paste0("**", space, "**"), space)) |>
# rename("sum loglik<sub>ICES rect</sub>" = rec,
# "sum loglik<sub>spatial</sub>" = space) |>
# dplyr::select(-rec_temp, -space_temp)
# AIC ices models
aic_ices <- AIC(fit_spr_m1, fit_her_m1, fit_sad_m1) |>
rownames_to_column() |>
mutate(Prey = "Sprat",
Prey = ifelse(rowname == "fit_her_m1", "Herring", Prey),
Prey = ifelse(rowname == "fit_sad_m1", "Saduria", Prey)) |>
mutate(AIC_r = AIC,
AIC_r2 = AIC) |>
dplyr::select(-rowname, -AIC, -df)
aic_space <- AIC(fit_spr_m2, fit_her_m2, fit_sad_m2) |>
rownames_to_column() |>
mutate(Prey = "Sprat",
Prey = ifelse(rowname == "fit_her_m2", "Herring", Prey),
Prey = ifelse(rowname == "fit_sad_m2", "Saduria", Prey)) |>
mutate(AIC_s = AIC,
AIC_s2 = AIC) |>
dplyr::select(-rowname, -AIC, -df)
aic <- aic_ices |>
left_join(aic_space, by = "Prey") |>
mutate(AIC_r = ifelse(AIC_r2 < AIC_s2, paste0("**", AIC_r, "**"), AIC_r),
AIC_s = ifelse(AIC_s2 < AIC_r2, paste0("**", AIC_s, "**"), AIC_s)) |>
rename("AIC<sub>spatial</sub>" = AIC_s,
"AIC<sub>ICES rect</sub>" = AIC_r) |>
dplyr::select(-AIC_s2, -AIC_r2)kableExtra::kbl(sum_loglik_space, format = "markdown", caption = "5-fold Spatial cross validation") |>
kable_styling(font_size = 20)| Prey | sum loglikICES rect | sum loglikspatial |
|---|---|---|
| Sprat | -17.230245 | -12.628093904989 |
| Herring | -8.682663 | -8.22049210072325 |
| Saduria | -15.292315 | -11.9475427099461 |
# kableExtra::kbl(sum_loglik_temp, format = "markdown", caption = "2-fold Temporal cross validation") |>
# kable_styling(font_size = 20)
kableExtra::kbl(aic, format = "markdown", caption = "AIC") |>
kable_styling(font_size = 20)| Prey | AICICES rect | AICspatial |
|---|---|---|
| Sprat | -4023.8311 | -4296.87560013517 |
| Herring | -265.7528 | -388.011502562052 |
| Saduria | -2128.6558 | -2271.95464531387 |
Check residuals and from selected models
# Summary
summary(fit_spr_m2)Spatial model fit by ML ['sdmTMB']
Formula: FR_spr ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
coef.est coef.se
depth_sc 0.33 0.07
spred_length_sc 0.66 0.08
Smooth terms:
Std. Dev.
sds(pred_length_sc) 32.42
Random intercepts:
Std. Dev.
month_f 0.65
Time-varying parameters:
coef.est coef.se
(Intercept)-1993 -6.19 0.42
(Intercept)-1994 -5.94 0.44
(Intercept)-1995 -7.15 0.40
(Intercept)-1996 -6.00 0.36
(Intercept)-1997 -7.00 0.43
(Intercept)-1998 -5.69 0.48
(Intercept)-1999 -6.08 0.35
(Intercept)-2000 -6.20 0.35
(Intercept)-2001 -5.35 0.37
(Intercept)-2002 -7.47 0.43
(Intercept)-2003 -4.62 0.37
(Intercept)-2004 -6.99 0.37
(Intercept)-2005 -5.83 0.35
(Intercept)-2006 -6.07 0.34
(Intercept)-2007 -6.23 0.34
(Intercept)-2008 -6.11 0.34
(Intercept)-2009 -5.77 0.35
(Intercept)-2012 -4.97 0.35
(Intercept)-2013 -5.57 0.34
(Intercept)-2014 -4.93 0.39
(Intercept)-2018 -6.22 0.42
(Intercept)-2021 -5.94 0.32
Dispersion parameter: 0.26
Tweedie p: 1.44
Matérn range: 19.81
Spatial SD: 1.20
ML criterion at convergence: -2157.438
See ?tidy.sdmTMB to extract these values as a data frame.
summary(fit_her_m2)Spatial model fit by ML ['sdmTMB']
Formula: FR_her ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
coef.est coef.se
depth_sc -0.13 0.09
spred_length_sc -0.43 0.07
Smooth terms:
Std. Dev.
sds(pred_length_sc) 38.1
Random intercepts:
Std. Dev.
month_f 0
Time-varying parameters:
coef.est coef.se
(Intercept)-1993 -8.97 0.50
(Intercept)-1994 -7.34 0.44
(Intercept)-1995 -6.75 0.36
(Intercept)-1996 -7.70 0.37
(Intercept)-1997 -5.99 0.34
(Intercept)-1998 -6.53 0.45
(Intercept)-1999 -6.92 0.31
(Intercept)-2000 -6.92 0.31
(Intercept)-2001 -6.82 0.33
(Intercept)-2002 -7.55 0.45
(Intercept)-2003 -5.51 0.29
(Intercept)-2004 -6.48 0.28
(Intercept)-2005 -7.45 0.30
(Intercept)-2006 -6.96 0.27
(Intercept)-2007 -7.22 0.28
(Intercept)-2008 -6.99 0.27
(Intercept)-2009 -6.39 0.30
(Intercept)-2012 -5.88 0.33
(Intercept)-2013 -5.74 0.28
(Intercept)-2014 -6.96 0.47
(Intercept)-2018 -7.66 0.42
(Intercept)-2021 -7.37 0.22
Dispersion parameter: 0.45
Tweedie p: 1.46
Matérn range: 17.38
Spatial SD: 1.41
ML criterion at convergence: -203.006
See ?tidy.sdmTMB to extract these values as a data frame.
**Possible issues detected! Check output of sanity().**
summary(fit_sad_m2)Spatial model fit by ML ['sdmTMB']
Formula: FR_sad ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
coef.est coef.se
depth_sc -0.44 0.15
spred_length_sc -0.04 0.09
Smooth terms:
Std. Dev.
sds(pred_length_sc) 12.28
Random intercepts:
Std. Dev.
month_f 0.64
Time-varying parameters:
coef.est coef.se
(Intercept)-1993 -8.64 0.93
(Intercept)-1994 -8.67 0.91
(Intercept)-1995 -8.37 0.87
(Intercept)-1996 -8.33 0.86
(Intercept)-1997 -8.05 0.86
(Intercept)-1998 -8.66 0.92
(Intercept)-1999 -8.34 0.84
(Intercept)-2000 -8.01 0.84
(Intercept)-2001 -7.30 0.84
(Intercept)-2002 -8.32 0.86
(Intercept)-2003 -7.84 0.84
(Intercept)-2004 -7.86 0.84
(Intercept)-2005 -7.38 0.83
(Intercept)-2006 -7.89 0.83
(Intercept)-2007 -8.54 0.84
(Intercept)-2008 -8.48 0.84
(Intercept)-2009 -8.37 0.85
(Intercept)-2012 -7.92 0.85
(Intercept)-2013 -8.24 0.84
(Intercept)-2014 -7.86 0.83
(Intercept)-2018 -7.44 0.79
(Intercept)-2021 -7.56 0.77
Dispersion parameter: 0.65
Tweedie p: 1.54
Matérn range: 83.70
Spatial SD: 2.59
ML criterion at convergence: -1144.977
See ?tidy.sdmTMB to extract these values as a data frame.
df |> filter(year == 2018) |> distinct(FR_sad)filter: removed 8,944 rows (97%), 265 rows remaining
distinct: removed 250 rows (94%), 15 rows remaining
# A tibble: 15 × 1
FR_sad
<dbl>
1 0.0368
2 0.0426
3 0.128
4 0
5 0.00219
6 0.00151
7 0.00400
8 0.00939
9 0.0218
10 0.0194
11 0.0740
12 0.00336
13 0.00233
14 0.00254
15 0.0115
# Residuals
spr_res <- mcmc_res <- residuals(fit_spr_m2, type = "mle-mcmc",
mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_spr_m2,
mcmc_iter = 201,
mcmc_warmup = 200))
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.002663 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 26.63 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 28.5231 seconds (Warm-up)
Chain 1: 0.05815 seconds (Sampling)
Chain 1: 28.5812 seconds (Total)
Chain 1:
df$spr_res <- as.vector(spr_res)
her_res <- mcmc_res <- residuals(fit_her_m2, type = "mle-mcmc",
mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_her_m2,
mcmc_iter = 201,
mcmc_warmup = 200))
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.002514 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 25.14 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 308.387 seconds (Warm-up)
Chain 1: 1.77732 seconds (Sampling)
Chain 1: 310.165 seconds (Total)
Chain 1:
Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
Warning: Examine the pairs() plot to diagnose sampling problems
df$her_res <- as.vector(her_res)
sad_res <- mcmc_res <- residuals(fit_sad_m2, type = "mle-mcmc",
mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_sad_m2,
mcmc_iter = 201,
mcmc_warmup = 200))
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.002502 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 25.02 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 66.2415 seconds (Warm-up)
Chain 1: 0.112717 seconds (Sampling)
Chain 1: 66.3542 seconds (Total)
Chain 1:
df$sad_res <- as.vector(sad_res)
# Plot all together
df |>
rename(Sprat = spr_res,
Herring = her_res,
Saduria = sad_res) |>
pivot_longer(c(Sprat, Herring, Saduria)) |>
ggplot(aes(sample = value)) +
facet_wrap(~name) +
stat_qq(size = 0.75, shape = 21, fill = NA) +
stat_qq_line() +
labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
theme(aspect.ratio = 1)rename: renamed 3 variables (Sprat, Herring, Saduria)
pivot_longer: reorganized (Sprat, Herring, Saduria) into (name, value) [was 9209x25, now 27627x24]
Plot conditional effects, random effects and spatial predictions
# Depth
vis_spr_dep <- visreg(fit_spr_m2, xvar = "depth_sc", plot = FALSE)
vis_her_dep <- visreg(fit_her_m2, xvar = "depth_sc", plot = FALSE)
vis_sad_dep <- visreg(fit_sad_m2, xvar = "depth_sc", plot = FALSE)
vis_dep <- bind_rows(vis_spr_dep$fit |> mutate(species = "Sprat"),
vis_her_dep$fit |> mutate(species = "Herring"),
vis_sad_dep$fit |> mutate(species = "Saduria")) |>
mutate(var = "Depth (scaled)") |>
rename(x = depth_sc)
d_dep <- bind_rows(vis_spr_dep$res |> mutate(species = "Sprat"),
vis_her_dep$res |> mutate(species = "Herring"),
vis_sad_dep$res |> mutate(species = "Saduria")) |>
mutate(var = "Depth (scaled)") |>
rename(x = depth_sc)
# Month
vis_spr_mon <- visreg(fit_spr_m2, xvar = "month_f", plot = FALSE)
vis_her_mon <- visreg(fit_her_m2, xvar = "month_f", plot = FALSE)
vis_sad_mon <- visreg(fit_sad_m2, xvar = "month_f", plot = FALSE)
vis_mon <- bind_rows(vis_spr_mon$fit |> mutate(species = "Sprat"),
vis_her_mon$fit |> mutate(species = "Herring"),
vis_sad_mon$fit |> mutate(species = "Saduria")) |>
mutate(var = "Month") |>
rename(x = month_f) |>
mutate(x = as.numeric(as.character(x)))
d_mon <- bind_rows(vis_spr_mon$res |> mutate(species = "Sprat"),
vis_her_mon$res |> mutate(species = "Herring"),
vis_sad_mon$res |> mutate(species = "Saduria")) |>
mutate(var = "Month") |>
rename(x = month_f) |>
mutate(x = as.numeric(as.character(x)))
# Predator length
vis_spr_len <- visreg(fit_spr_m2, xvar = "pred_length_sc", plot = FALSE)
vis_her_len <- visreg(fit_her_m2, xvar = "pred_length_sc", plot = FALSE)
vis_sad_len <- visreg(fit_sad_m2, xvar = "pred_length_sc", plot = FALSE)
vis_len <- bind_rows(vis_spr_len$fit |> mutate(species = "Sprat"),
vis_her_len$fit |> mutate(species = "Herring"),
vis_sad_len$fit |> mutate(species = "Saduria")) |>
mutate(var = "Predator length") |>
rename(x = pred_length_sc)
d_len <- bind_rows(vis_spr_len$res |> mutate(species = "Sprat"),
vis_her_len$res |> mutate(species = "Herring"),
vis_sad_len$res |> mutate(species = "Saduria")) |>
mutate(var = "Predator length") |>
rename(x = pred_length_sc)
vis <- bind_rows(vis_dep, vis_mon, vis_len)
vis_dat <- bind_rows(d_dep, d_mon, d_len)
# Plot!
ggplot(vis |> filter(!var == "Month"), aes(x = x, y = visregFit)) +
facet_grid(var ~ species) +
geom_point(data = vis_dat |> filter(!var == "Month"), aes(x = x, y = visregRes),
alpha = 0.3, color = "gray50") +
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.3, color = NA) +
geom_line(color = "steelblue", linewidth = 1) +
labs(x = "Scaled variable", y = "Prediction") +
NULLggsave(paste0(home, "/figures/conditional.pdf"), width = 17, height = 11, units = "cm")
# Now do month (categorial)
ggplot(vis |> filter(var == "Month"), aes(x = as.factor(x), y = visregFit)) +
facet_wrap(~ species) + # free grid?
geom_jitter(data = vis_dat |> filter(var == "Month"), aes(x = as.factor(x), y = visregRes),
alpha = 0.1, color = "gray70") +
geom_errorbar(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.6, color = "steelblue",
width = 0, linewidth = 1) +
geom_point(color = "steelblue", size = 2) +
labs(x = "Month", y = "Prediction") +
theme(aspect.ratio = 5/6) +
NULLggsave(paste0(home, "/figures/conditional_month.pdf"), width = 17, height = 6, units = "cm")Calculate indices
# Need to refit the models with the extra time argument. The reason we don't do that before is because it creates a mismatch in residual dimensions and data
# This is for interpolating between year using the random walk
extra_time <- pred_grid |> filter(!year %in% df$year) |> distinct(year) |> pull(year)filter: removed 357,852 rows (76%), 113,862 rows remaining
distinct: removed 113,855 rows (>99%), 7 rows remaining
# Spatial random field
fit_spr_m2_extra <- sdmTMB(
FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
extra_time = extra_time,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))
fit_her_m2_extra <- sdmTMB(
FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
extra_time = extra_time,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))
fit_sad_m2_extra <- sdmTMB(
FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
extra_time = extra_time,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))
# Predict on grid, for indices and maps
pred_spr <- predict(fit_spr_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)
pred_her <- predict(fit_her_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)
pred_sad <- predict(fit_sad_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)
# Make temporal index!
ncells <- filter(pred_grid, year == max(pred_grid$year)) |> nrow()filter: removed 455,448 rows (97%), 16,266 rows remaining
index_spr <- get_index(pred_spr, area = rep(1/ncells, nrow(pred_spr$data)))Bias correction is turned off.
It is recommended to turn this on for final inference.
index_her <- get_index(pred_her, area = rep(1/ncells, nrow(pred_her$data)))Bias correction is turned off.
It is recommended to turn this on for final inference.
index_sad <- get_index(pred_sad, area = rep(1/ncells, nrow(pred_sad$data)))Bias correction is turned off.
It is recommended to turn this on for final inference.
# Make long
index <- bind_rows(index_spr |> mutate(prey = "Sprat"),
index_her |> mutate(prey = "Herring"),
index_sad |> mutate(prey = "Saduria")) |>
mutate(Observed = ifelse(year %in% extra_time, "No", "Yes"))mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'Observed' (character) with 2 unique values and 0% NA
# For comparison with data
df_sum <- df |>
group_by(year) |>
summarise(Sprat = mean(FR_spr),
Herring = mean(FR_her),
Saduria = mean(FR_sad)) |>
pivot_longer(c("Sprat", "Herring", "Saduria"), names_to = "prey", values_to = "mean_fr")group_by: one grouping variable (year)
summarise: now 22 rows and 4 columns, ungrouped
pivot_longer: reorganized (Sprat, Herring, Saduria) into (prey, mean_fr) [was 22x4, now 66x3]
sort(unique(index$year)) [1] 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
[16] 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021
sort(unique(pred_grid$year)) [1] 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
[16] 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021
index |>
ggplot(aes(year, est, fill = Observed)) +
geom_point(shape = 21, alpha = 0.7) +
scale_fill_manual(values = c("white", "grey10")) +
facet_wrap(~prey, ncol = 1, scales = "free") +
geom_errorbar(aes(ymin = lwr, ymax = upr), alpha = 0.4, width = 0) +
geom_point(data = df_sum, aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
labs(x = "Year", y = "Feeding ratio", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 3)) +
theme(legend.position = "bottom")ggsave(paste0(home, "/figures/supp/index_ci.pdf"), width = 11, height = 20, units = "cm")
index |>
mutate(upr = ifelse(prey == "Sprat" & upr > 0.015, 0.015, upr)) |>
ggplot(aes(year, est)) +
geom_point(alpha = 0.7) +
stat_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue") +
facet_wrap(~prey, ncol = 3, scales = "free") +
labs(x = "Year", y = "Feeding ratio", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.9, 0.9))mutate: changed 10 values (11%) of 'upr' (0 new NA)
ggsave(paste0(home, "/figures/index.pdf"), width = 17, height = 6, units = "cm")Plot spatial predictions
spatial_preds <- bind_rows(pred_spr$data |> mutate(prey = "Sprat"),
pred_her$data |> mutate(prey = "Herring"),
pred_sad$data |> mutate(prey = "Saduria"))mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
# Plot spatial random effect
plot_map +
geom_raster(data = spatial_preds |> filter(year == 2000),
aes(X*1000, Y*1000, fill = omega_s)) +
scale_fill_gradient2(name = "Spatial random field") +
facet_wrap(~prey) +
geom_sf(size = 0.1) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.05, 0.65)) +
theme(legend.position = c(0.05, 0.7),
legend.key.height = unit(0.6, "line"),
legend.key.width = unit(0.2, "line"))filter: removed 1,366,344 rows (97%), 48,798 rows remaining
Warning: Removed 2283 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/diet_omega_s.pdf"), width = 17, height = 6, units = "cm")Warning: Removed 2283 rows containing missing values (`geom_raster()`).
# Plot spatiotemporal predictions
plot_map_fc +
geom_raster(data = spatial_preds |> filter(prey == "Sprat"), aes(X*1000, Y*1000, fill = exp(est))) +
scale_fill_viridis(trans = "sqrt", name = "FR") +
facet_wrap(~year) +
geom_sf(size = 0.1)filter: removed 943,428 rows (67%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/sprat_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")Warning: Removed 22069 rows containing missing values (`geom_raster()`).
plot_map_fc +
geom_raster(data = spatial_preds |> filter(prey == "Herring"), aes(X*1000, Y*1000, fill = exp(est))) +
scale_fill_viridis(trans = "sqrt", name = "FR") +
facet_wrap(~year) +
geom_sf(size = 0.1)filter: removed 943,428 rows (67%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/herring_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")Warning: Removed 22069 rows containing missing values (`geom_raster()`).
plot_map_fc +
geom_raster(data = spatial_preds |> filter(prey == "Saduria"), aes(X*1000, Y*1000, fill = exp(est))) +
scale_fill_viridis(trans = "sqrt", name = "FR") +
facet_wrap(~year) +
geom_sf(size = 0.1)filter: removed 943,428 rows (67%), 471,714 rows remaining
Warning: Removed 22069 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/saduria_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")Warning: Removed 22069 rows containing missing values (`geom_raster()`).
# Plot spatiotemporal predictions for a year and all species
plot_map +
geom_raster(data = spatial_preds|> filter(year == "2000"), aes(X*1000, Y*1000, fill = exp(est))) +
scale_fill_viridis(trans = "sqrt", name = "FR") +
facet_wrap(~prey) +
geom_sf(size = 0.1) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.05, 0.7),
legend.key.height = unit(0.6, "line"),
legend.key.width = unit(0.2, "line"))filter: removed 1,366,344 rows (97%), 48,798 rows remaining
Warning: Removed 2283 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/spatial_prey_prediction_2000.pdf"), width = 17, height = 6, units = "cm")Warning: Removed 2283 rows containing missing values (`geom_raster()`).
Plot correlations between spatial overlap and population-level FR
# Join feeding ratio indicies (per capita) and overlap
cod_spr <- read_csv(paste0(home, "/output/cod_pel_sum_ovrlap.csv")) |>
dplyr::select(year, cod_spr_ovr_tot) |>
rename(overlap = cod_spr_ovr_tot) |>
mutate(prey = "Sprat")Rows: 27 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): year, cod_spr_ovr_tot, cod_spr_ovr_sd_tot, cod_her_ovr_tot, cod_her...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)
mutate: new variable 'prey' (character) with one unique value and 0% NA
cod_her <- read_csv(paste0(home, "/output/cod_pel_sum_ovrlap.csv")) |>
dplyr::select(year, cod_her_ovr_tot) |>
rename(overlap = cod_her_ovr_tot) |>
mutate(prey = "Herring")Rows: 27 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): year, cod_spr_ovr_tot, cod_spr_ovr_sd_tot, cod_her_ovr_tot, cod_her...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)
mutate: new variable 'prey' (character) with one unique value and 0% NA
cod_ben <- read_csv(paste0(home, "/output/cod_ben_sum_ovrlap.csv")) |>
dplyr::select(year, cod_sad_ovr_tot) |>
rename(overlap = cod_sad_ovr_tot) |>
mutate(prey = "Saduria")Rows: 54 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (3): year, quarter, cod_sad_ovr_tot
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)
mutate: new variable 'prey' (character) with one unique value and 0% NA
overlap <- bind_rows(cod_spr, cod_her, cod_ben)
index_ovr <- index |>
left_join(overlap, by = c("year", "prey")) |>
drop_na()left_join: added one column (overlap)
> rows only in x 6
> rows only in y ( 0)
> matched rows 108 (includes duplicates)
> =====
> rows total 114
drop_na: removed 6 rows (5%), 108 rows remaining
# Join in cod biomass data to calculate snapshot predation
cod_pred <- read_csv(paste0(home, "/output/pred_cod.csv")) |>
filter(quarter == 4)Rows: 943428 Columns: 38
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): substrate, month_year, ices_rect
dbl (35): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_d...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 471,714 rows (50%), 471,714 rows remaining
# Left join cod biomass density onto the grid of diet predictions
spatial_preds_dens <- left_join(spatial_preds, cod_pred |>
dplyr::select(est, X, Y, year) |>
rename(est_codbiom = est))rename: renamed one variable (est_codbiom)
Joining with `by = join_by(X, Y, year)`left_join: added one column (est_codbiom)
> rows only in x 0
> rows only in y ( 0)
> matched rows 1,415,142
> ===========
> rows total 1,415,142
# Multiply local feeding ration with biomass density of cod
spatial_preds_dens <- spatial_preds_dens |>
mutate(pred = exp(est) * est_codbiom)mutate: new variable 'pred' (double) with 1,415,142 unique values and 0% NA
# Summarise across years
spatial_preds_dens_sum <- spatial_preds_dens |>
group_by(prey, year) |>
summarise(tot_pred = sum(pred))group_by: 2 grouping variables (prey, year)
summarise: now 87 rows and 3 columns, one group variable remaining (prey)
# Left_join with index data
index_ovr <- index_ovr |>
left_join(spatial_preds_dens_sum, by = c("year", "prey")) |>
drop_na()left_join: added one column (tot_pred)
> rows only in x 0
> rows only in y ( 6)
> matched rows 108
> =====
> rows total 108
drop_na: no rows removed
# Check some regressions...
summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Sprat")))filter: removed 81 rows (75%), 27 rows remaining
Call:
lm(formula = tot_pred ~ overlap, data = filter(index_ovr, prey ==
"Sprat"))
Residuals:
Min 1Q Median 3Q Max
-21096 -8859 -3394 4497 50120
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3876 13500 -0.287 0.776
overlap 59012 41360 1.427 0.166
Residual standard error: 15530 on 25 degrees of freedom
Multiple R-squared: 0.0753, Adjusted R-squared: 0.03831
F-statistic: 2.036 on 1 and 25 DF, p-value: 0.166
summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Herring")))filter: removed 81 rows (75%), 27 rows remaining
Call:
lm(formula = tot_pred ~ overlap, data = filter(index_ovr, prey ==
"Herring"))
Residuals:
Min 1Q Median 3Q Max
-4959.6 -3374.5 -2305.8 468.6 13247.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3900 4986 0.782 0.441
overlap 3804 9792 0.388 0.701
Residual standard error: 5332 on 25 degrees of freedom
Multiple R-squared: 0.006, Adjusted R-squared: -0.03376
F-statistic: 0.1509 on 1 and 25 DF, p-value: 0.701
summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Saduria")))filter: removed 54 rows (50%), 54 rows remaining
Call:
lm(formula = tot_pred ~ overlap, data = filter(index_ovr, prey ==
"Saduria"))
Residuals:
Min 1Q Median 3Q Max
-4230.7 -2127.8 -73.2 1827.6 7535.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1140 1388 0.821 0.415
overlap 12247 3512 3.487 0.001 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2719 on 52 degrees of freedom
Multiple R-squared: 0.1896, Adjusted R-squared: 0.174
F-statistic: 12.16 on 1 and 52 DF, p-value: 0.001001
summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Sprat")))filter: removed 81 rows (75%), 27 rows remaining
Call:
lm(formula = est ~ overlap, data = filter(index_ovr, prey ==
"Sprat"))
Residuals:
Min 1Q Median 3Q Max
-0.0049434 -0.0014152 -0.0004706 0.0008664 0.0080607
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0002763 0.0024015 -0.115 0.9093
overlap 0.0145854 0.0073571 1.983 0.0585 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.002762 on 25 degrees of freedom
Multiple R-squared: 0.1359, Adjusted R-squared: 0.1013
F-statistic: 3.93 on 1 and 25 DF, p-value: 0.05851
summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Herring")))filter: removed 81 rows (75%), 27 rows remaining
Call:
lm(formula = est ~ overlap, data = filter(index_ovr, prey ==
"Herring"))
Residuals:
Min 1Q Median 3Q Max
-0.0008822 -0.0005791 -0.0003672 0.0003922 0.0025145
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0005373 0.0008206 0.655 0.519
overlap 0.0012624 0.0016116 0.783 0.441
Residual standard error: 0.0008775 on 25 degrees of freedom
Multiple R-squared: 0.02396, Adjusted R-squared: -0.01509
F-statistic: 0.6136 on 1 and 25 DF, p-value: 0.4408
summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Saduria")))filter: removed 54 rows (50%), 54 rows remaining
Call:
lm(formula = est ~ overlap, data = filter(index_ovr, prey ==
"Saduria"))
Residuals:
Min 1Q Median 3Q Max
-0.0007962 -0.0005396 -0.0002388 0.0002035 0.0019474
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0014897 0.0003855 3.865 0.00031 ***
overlap -0.0003224 0.0009749 -0.331 0.74222
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.0007548 on 52 degrees of freedom
Multiple R-squared: 0.002098, Adjusted R-squared: -0.01709
F-statistic: 0.1093 on 1 and 52 DF, p-value: 0.7422
# Calculate correlations between FR, predation and overlap
cor <- plyr::ddply(index_ovr, c("prey"),
summarise,
cor_fr_ovr = round(cor(overlap, est), 2),
cor_pred_ovr = round(cor(overlap, tot_pred), 2)) |>
pivot_longer(c("cor_fr_ovr", "cor_pred_ovr")) |>
mutate(name = ifelse(name == "cor_fr_ovr", "Feeding ratio", "Total predation"))summarise: now one row and 2 columns, ungrouped
summarise: now one row and 2 columns, ungrouped
summarise: now one row and 2 columns, ungrouped
pivot_longer: reorganized (cor_fr_ovr, cor_pred_ovr) into (name, value) [was 3x3, now 6x3]
mutate: changed 6 values (100%) of 'name' (0 new NA)
# Plot!
library(ggh4x)<
index_ovr |>
rename('Feeding ratio' = est,
'Total predation' = tot_pred) |>
pivot_longer(c('Feeding ratio', 'Total predation')) |>
ggplot(aes(overlap, value)) +
geom_point() +
#facet_wrap(name ~ prey, scales = "free") +
ggh4x::facet_grid2(name ~ prey, scales = "free", independent = "y") +
geom_smooth(aes(overlap, value), inherit.aes = FALSE, color = "steelblue",
formula = y~s(x, k=10), method = "gam") +
geom_text(data = cor,
aes(label = paste("r=", value, sep = "")), x = -Inf, y = Inf, hjust = -.1, vjust = 2.5,
inherit.aes = FALSE, fontface = "italic", size = 2.5, color = "tomato3") +
labs(y = "Population-level feeding ratio", x = "Spatial overlap", color = "Year") +
theme_sleek(base_size = 9) +
theme(legend.position = "bottom", aspect.ratio = 1) +
NULLrename: renamed 2 variables (Feeding ratio, Total predation)
pivot_longer: reorganized (Feeding ratio, Total predation) into (name, value) [was 108x10, now 216x10]
Warning in library(ggh4x) < ggplot(pivot_longer(rename(index_ovr, `Feeding
ratio` = est, : longer object length is not a multiple of shorter object length
[1] TRUE TRUE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE TRUE FALSE FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE
[25] FALSE FALSE FALSE TRUE TRUE FALSE TRUE FALSE FALSE
# Summarise the differences here, make fake facet grid free axis
ggsave(paste0(home, "/figures/fr_overlap_cor.pdf"), width = 17, height = 11, units = "cm")TODO:
- there are new data in the db that are not in stefans data also for the older period. find a way to not get duplicates.
- update as I get more data
- explore the spatial effects more
- make a clean, straight script for all species that produces the plots for the report